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ABSTRACT 

Context. The nearby Large Magellanic Cloud (LMC) provides a rare opportunity of a spatially resolved view of an external star-forming galaxy 
in y-rays. The LMC was detected at 0.1-100 GeV as an extended source with CGRO/EGRET and using early observations with the Ferntz-LAT. 
The emission was found to correlate with massive star-forming regions and to be particularly bright towards 30 Doradus. 

Aims. Studies of the origin and transport of cosmic rays (CRs) in the Milky Way are frequently hampered by line-of-sight confusion and poor 
distance determination. The LMC offers a complementary way to address these questions by revealing whether and how the y-ray emission is 
connected to specific objects, populations of objects, and structures in the galaxy. 

Methods. We revisited the y-ray emission from the LMC using about 73 months of Fermi-hAT P7REP data in the 0.2-100 GeV range. We devel¬ 
oped a complete spatial and spectral model of the LMC emission, for which we tested several approaches: a simple geometrical description, 
template-fitting, and a physically driven model for CR-induced interstellar emission. 

Results. In addition to identifying PSR J0540-6919 through its pulsations, we find two hard sources positionally coincident with plerion N 157B 
and supernova remnant N 132D, which were also detected at TeV energies with H.E.S.S. We detect an additional soft source that is currently 
unidentified. Extended emission dominates the total flux from the LMC. It consists of an extended component of about the size of the galaxy and 
additional emission from three to four regions with degree-scale sizes. If it is interpreted as CRs interacting with interstellar gas, the large-scale 
emission implies a large-scale population of ~ 1-100 GeV CRs with a density of ~30% of the local Galactic value. On top of that, the three to four 
small-scale emission regions would correspond to enhancements of the CR density by factors 2 to 6 or higher, possibly more energetic and younger 
populations of CRs compared to the large-scale population. An alternative explanation is that this is emission from an unresolved population of 
at least two dozen objects, such as pulsars and their nebulae or supernova remnants. This small-scale extended emission has a spatial distribution 
that does not clearly correlate with known components of the LMC, except for a possible relation to cavities and supergiant shells. 

Conclusions. The Fernii-hKI GeV observations allowed us to detect individual sources in the LMC. Three of the newly discovered sources are 
associated with rare and extreme objects. The 30 Doradus region is prominent in GeV y-rays because PSR J0540-6919 and N 157B are strong 
emitters. The extended emission from the galaxy has an unexpected spatial distribution, and observations at higher energies and in radio may help 
to clarify its origin. 

Key words. Gamma rays: galaxies - Galaxies: individual: Large Magellanic Cloud - Cosmic rays 


1. Introduction 

The Large Magellanic Cloud (LMC) is a unique target for high- 
energy astrophysics. Given the sensitivity and angular resolution 
of current y-ray instruments, it is the best opportunity we have to 
obtain an external and spatially resolved view of a star-forming 


galaxy. The Small Magellanic Cloud (SMC) offers a similar op¬ 
portunity, but the galaxy is smaller, farther away, less active 
in terms of star for mation, and has a mor e complex and less 
favourable geometry (IScowcroft et alJ2015l) . Other star-forming 
galaxies within reach of current y-ray instruments are at least an 
order of magnitude more distant (e.g. M31, M82, and NGC 253; 
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lAbdo et ^l2010all3lAbramowski et'^l2012h and have not been 
spatially resolved (in either the GeV or TeV domain). 

In the Milky Way, hundreds of discrete objects of differ¬ 
ent classes have been detected, together with diffuse interstellar 
emission over a broad range of spatial scales. Studies of Galactic 
sources are often complicated by line-of-sight confusion and in¬ 
accurate distance estimates, however. In contrast, the LMC of¬ 
fers more modest prospects in terms of the number of detectable 
sources, but it provides an outside and global perspective on a 
whole galaxy seen at low inclination, whose dist ance is deter¬ 
mined at the percent level (iPietrzvnski et al.lF2013h . 

Studying the LMC in y-rays can be a valuable complement 
to Milky Way studies on a variety of topics: particle accelerators 
such as supernova remnants (SNRs), pulsars, and pulsar wind 
nebulae (PWNe), with the possibility to find more rare and ex¬ 
treme objects; diffuse interstellar emission, which may show us 
how a cosmic-ray population developed in another galaxy and 
thus may be a useful test of our current theories on cosmic-ray 
transport; indirect searches for dark matter, which provides lim¬ 
its competitive with recent c onstraints, and that h ave different 
uncertainties and systematics (iBucklev et al.ll201~5h . On the first 
two subjects, the outside point of view may reveal how cosmic- 
ray sources and transport are related to the structure and content 
of a galaxy as a whole. Here, building a global picture will be 
facilitated by the fact that the LMC was deeply surveyed over 
the past decades by almost all instruments that were able to do 
so, from radio to y-rays. 

From 11 months of Fermi-LA J' continuous all sky-survey 
observations, lAbdo et ^ (1201 Od) . hereafter Paper I, reported 
the 33cr detection of the LMC in y-rays. The source was spa¬ 
tially extended and its spectrum was consistent with the ob¬ 
served y-rays originating from cosmic rays (CRs) interacting 
with the interstellar medium (ISM) through inverse-Compton 
scattering, Bremsstrahlung, and hadronic interactions, but con¬ 
tributions from discrete objects such as pulsars were not ruled 
out at that time. The emission was found to be relatively strong 
in the direction of the 30 Doradus star-forming region; more gen¬ 
erally, the emission seemed spatially correlated with tracers of 
massive star-forming regions, such as the Her emission or the 
population of known Wolf-Rayet stars, and showed little corre¬ 
lation with the gas column density distribution. This was inter¬ 
preted as evidence supporting the idea that CRs are accelerated 
in massive star-forming regions as a result of the large amounts 
of kinetic energy released by the stellar winds and supernova 
explosions of massive stars. In addition, the close confinement 
of y-ray emission to star-forming regions suggested a relatively 
short GeV CR proton diffusion length. 

We now have over six times more data than were used in 
iPaner ll and better instrument performance. In this paper, we 
present a new analysis of the GeV y-ray emission from the LMC. 
The paper is organised as follows: we first introduce the data set 
used in this work (Sect. El and the models developed to account 
for the emission seen in the LMC region (Sect.Ei then, the prop¬ 
erties of the point-like sources found within the LMC boundaries 
are introduced and their identification or association with known 
objects is discussed (Sect. E; last, the extended y-ray emission 
from the LMC is presented, and a possible interpretation in terms 
of CR-induced interstellar emission or unresolved population of 
sources is given (Sect.E- Throughout t he paper, the distance o f 
the LMC is assumed to be c/ = 50.1 kpc (IPietrzvnski et al.l2013h . 
and its incl ination with respect to the line of sight is assumed to 
be i — 30° (Ivan der Marelll 206 ^ . 


2. Data selection 

The following analysis was performed using 73.3 months 
of observations with the Fermi LAT (mission elapsed time 
239587200 to 43269 4964) , primarily taken in all-sky survey 
mode (lAtwood et al.l 1200^ . The data set was produced with 
the so-called Pass 7 reprocessed (P7REP) version of the event 
analysis and selection criteria, which takes into account effects 
measured in flight that were not considered in pre-launch per¬ 
formance estimates, such as pile- up and accidental coincid ence 
effects in the detector subsystems (lAckermann et al.ll2012al) . and 
updated calibration constants, to includ e effects such as th e 
degradation in the calorimeter light yield (iBregeon et al.ll2()T^ . 
SOURCE class events were used in this work, excluding those 
coming from zenith angles larger than 100° or detected when the 
rocking angle of the satellite was larger than 52° to reduce con¬ 
tamination by atmospheric y-rays from the Earth. We analysed 
together events converting in the front and back sections of the 
LAT, but we checked that considering them separately does not 
alter the results. The analysis was performed using Science Tools 
package v09r32p05 and using the P7REP_SOURCE_V15 instru¬ 
ment response functions (IRFs), all available from the Fermi 
Science Support Centefl 

A region of interest (ROI) specific to the LMC was defined 
as a 10° X 10° square centred on (a,S) - (80.894°,-69.756°) 
and aligned on equatorial coordinates (here and throughout the 
paper, equatorial coordinates correspond to the J2000.0 epochfl 
Events contained in this region were selected from the data set 
described above. The energy range considered in this analysis 
is 0.2-100 GeV; the lower energy bound is dictated by the poor 
angular resolution at the lowest energies. A counts map of the 
ROI is shown in the top panel of Fig. E 

3. Emission modelling 

We aim at building a complete model for the spatial and spectral 
distribution of the counts in the ROI and energy range defined 
above. As summarised below, we tested several possible ap¬ 
proaches, each providing different insights into the y-ray emis¬ 
sion properties. All rely on a model-fitting procedure following 
a maximum likelihood approach for binned data and Poisson 
statistics. A given model consists of several emission compo¬ 
nents and has a certain number of free parameters. A distribu¬ 
tion of expected counts in position and energy is obtained by 
convolving the model with the IRFs, taking into account the ex¬ 
posure achieved for the data set that is used. The free param¬ 
eters are then adjusted in an iterative way until the distribu¬ 
tion of expected counts provides the highest likelihood of the 
data given the model. Fits were made over the 0.2-100 GeV 
range with 17 logarithmic energy bins and using 0.1° x 0.1° pix¬ 
els. In the following, models are compared based on the maxi¬ 
mum value of the logarithm of the likelihood function, denoted 
log X- The significance of model components or additional pa¬ 
rameters is evaluated using the test statistic, whose expression is 
TS - 2(log X - log Xo), where Lq is the likelihood of the ref¬ 
erence model without the additional parameter or component. A 
summary of the likelihoods and numbers of degrees of freedom 
of each model is given in Sect. 13.51 


* See http://fermi.gsfc.nasa.gov/ssc/ 

^ As reported below, using a larger ROI had no impact on the results. 
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Fig. 1. Total counts map in the 0.2-100 GeV band and residual 
counts map after subtracting the background model described in 
Sect. B.lb top and bottom panels, respectively). Both maps have 
0.1° X 0.1° pixels and were smoothed with a Gaussian kernel 
with cr = 0.2°. Colours are displayed on a square-root scale. 
White lines are contours of the atomic hydrogen distribution in 
the LMC at a relative value of 1/8 of the peak in the distribution 
(see Sect. l3.3l for the origin of the data). 


3.1. Background model 

As a first step in the process of modelling the emission over the 
ROI and before developing a model for the LMC, we have to 
account for known background and foreground emission in the 
form of diffuse or isolated sources. These are: 

1- The Galactic interstellar emission, arising from CRs in¬ 
teracting with the ISM in our Galaxy. In the Fenni-LAI energy 
range, this e mission is dominated by h adronic emission from in¬ 
terstellar gas (lAckermann et al.ll2012ch . Even at the Galactic lat¬ 


itude of the LMC, b ~ -33°, this foreground radiation is clearly 
present as structured emission in the counts map. We modelled 
it using the template provided by the Fermi Science Support 
Center, gllJem_v05_revl.fit, with free normalisatio n in the fit. 
In the preparation of this template (ICasandiia BUMl, any signal 
from gas in the LMC was removed, so that the y-ray emission 
from the LMC is not erroneously absorbed in the Galactic diffuse 
emission model. Importantly, the LMC region is not affected by 
the large-scale residual structures re-injected into the template. 

2- An isotropic background, which accounts for an approxi¬ 
mately isotropic diffuse y-ray emission component and residual 
CRs misclassified as y-rays in the LAT. The origin of the as- 
trophysical emission is currently unclear and it may come from 
multiple s ources, ranging from the s olar system to cosmological 
structures (lAckermann et al.ll2015al) . It was modelled using the 
publicly available isotropic spectral model isojource-v05.txt, 
with free normalisation in the fit. 

3- All isolated sources in the region that were previously 
detected and listed in the F enni-LP^ second source catalogue 
(2FGL, iNolan et'aDl2012h . A total of seven sources fall into 
the ROI defined above. Four of them were dismissed because 
they are located within the LMC boundaries and may ac¬ 
tually correspond to components of the LMC emission that 
we aimed at modelling. Source 2FGL J0532.5-7223 was ex¬ 
cluded because its significance was too low (below 3cr), but 
a source not listed in the 2FGL catalogue was found nearby, 
at the position (a, 6) - (82.4°,-72.7°), with a TS above 100, 
a power-law spectrum, and a variable flux. This additional 
source is present in the Fermi-LAiY th ird source catalogue as 
3FGL J0529.8-7242 (lAcero et al.l2015l) . The other field sources 
are 2FGL J0438.0-7331 and 2FGL J0601.1 -7037, the latter be¬ 
ing associated with the radio source PKS 0601-70 and exhibit¬ 
ing strongly variable emission. All three were included in the 
model as point-like sources using the spectral shapes identified 
as most suitable in the 2FGL catalogue and leaving their spectral 
parameters free in the fit. We also included sources lying outside 
the ROI, up to a distance of 3°, to account for spillover of their 
emission inside the ROI at low energies, where the point-spread 
function has a degree-scale size. A total of ten such sources were 
included in the model, with spectral shapes and parameters fixed 
at the catalogue values. 

All the components described above form the basis of the 
emission model and are referred to as the background model. 
It has a total of nine degrees of freedom (one for each diffuse 
emission template, three for 2FGL J060L1-7037, and two each 
for 2FGL J0438.0-7331 and 3FGL J0529.8-7242). We now de¬ 
scribe the modelling of the excess signal that is not accounted 
for by this background model. 

3.2. Analytic model 

Starting from the background model, we first aim to describe 
the remaining emission with a combination of point-like and 
2D Gaussian-shaped spatial intensity distributions, adding new 
components successively. 

Point-like sources can be identified if they have hard spectra 
and are bright enough, because the angular resolution at high en¬ 
ergies >10-20 GeV is relatively good and allows distinguishing 
them from any extended emission. Inspection of the >20 GeV 
counts map suggested the presence of two such sources (called 
P2 and P4 in the following), the significance and point-like na¬ 
ture of which was confirmed by subsequent analyses. In addi¬ 
tion, one source was identified as a y-ray pulsar from its char¬ 
acteristic pulsations: the source called PI in the following was 
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Table 1. Point sources found in the LMC 


Name 

Coordinates 

Error radius 

Spectrum 

TS 

PI 

85.0465°,-69.3316° 


PEC 

151 

P2 

84.43°, -69.17° 

0.02° 

PL 

96 

P3 

83.54°, -67.54° 

0.06° 

PL 

116 

P4 

81.13°,-69.62° 

0.03° 

PL 

28 


From left to right, the columns list the name of the point source (PI being firmly identified as PSR J0540-6919), equatorial coordinates, error 
on the position, best-fit spectral shape (PL=power law, PLC=power law with exponential cutoff), and the test statistic of the source. All values 
correspond to the emissivity model. Position uncertainties correspond to a 95% confidence level. 


Table 2. Emission components of the analytic model 


Name 

Coordinates 

Error radius 

Size 

Spectrum 

TS 

G1 

79.70°, -68.55° 

0.2° 

1.95° ±0.05° 

LP 

905 

G2 

83.40°,-69.15° 

0.1° 

0.70° ± 0.05° 

PL 

251 

G3 

83.10°,-66.60° 

0.2° 

0.45° ±0.10° 

PL 

211 

G4 

74.20°, -69.50° 

0 

d 

0.40° ±0.10° 

PL 

55 


From left to right, the columns list the name of the 2D Gaussian intensity distribution component, equatorial coordinates of the centre of the 2D 
Gaussian, error on the position, cr parameter of the 2D Gaussian, best-fit spectral shape (PL=power law, LP=log-parabola), and the test statistic of 
the emission component. Position uncertainties correspond to a 95% confidence level, size uncertainties correspond to Icr. 


Table 3. Emissivity components of the emissivity model 


Name 

Coordinates 

Error radius 

Size 

Spectrum 

TS 

EO 

80.00°, -68.00° 

0.5° 

d 

0 

LP 

891 

El 

82.40°, -68.85° 

0.4° 

0.6° ±0.1° 

PL 

276 

E2 

82.97°, -66.65° 

0.2° 

0.4° ±0.1° 

PL 

270 

E3 

82.25°, -69.25° 

0.2° 

0.3° ±0.1° 

PL 

276 

E4 

75.25°, -69.75° 

0.4° 

0.6° ± 0.2° 

PL 

120 


From left to right, the columns list the name of the emissivity component, equatorial coordinates of the centre of the 2D Gaussian emissivity dis¬ 
tribution, error on the position, <r parameter of the 2D Gaussian emissivity distribution, best-fit spectral shape (PL=power law, LP=log-parabola), 
and the test statistic of the emission component (the TS for El and E3 is that obtained when both components are tied together, see text). Position 
uncertainties correspond to a 95% confidence level, size uncertainties correspond to Icr. 


identified as PSR J0540-6919 (see Sect. I4.l| i. and its position 
was fixed to the posit ion of the pulsar from optical observations 
(iMignani et al.ll20lO) . Starting from these three point-like con¬ 
tributions to the LMC emission, an iterative procedure was used 
to characterise the remaining emission. 

At each step, a scan over position and size of the new emis¬ 
sion component is performed to identify the source that provides 
the best fit to the data. Eor each trial position and size, the full 
model is multiplied by the exposure, convolved with the point 
spread function, and fit to the data in a binned maximum likeli¬ 
hood analysis for Poisson statistics. A power-law spectral shape 
is assumed for the new component (which turned out to be a 
good approximation for most components). If the improvement 
of the likelihood is significant - which we defined as TS >25 
- the component is added to the model and a new iteration 
starts. The process stops when adding a new component yields a 
TS <25. The subsequent step is to re-optimize the positions and 
sizes of the components, from the brightest to the faintest in turn. 
The final stage is deriving bin-by-bin spectra for all components 
to check that the initially adopted power-law spectral shape is 
appropriate. If not, it is replaced by a power-law with exponen¬ 
tial cutoff or a log-parabola shape, depending on which provides 
the best fit and a significant improvement over the power-law 
ass umption (for the fo rmulae of the different spectral functions, 
see iNolan etalJl2012h . 


This procedure resulted in an emission model with eight 
components: four point-like objects and four Gaussian-shaped, 
spatially extended components. In the following, this model is 
referred to as the analytic model. The properties of all compo¬ 
nents are summarised in Tables[T]and|3 and their layout is illus¬ 
trated in Eig.|2](the positions listed in Table [1] correspond to the 
emissivity model described below; the positions corresponding 
to the analytic model are very similar and fully consistent with 
them, as illustrated in Eig. H]). Point sources are labelled PI to 
P4, and extended sources are labelled G1 to G4. 

3.3. Emissivity model 

The diffuse emission part of the analytic model includes a large- 
scale contribution, which is located close to the geometrical cen¬ 
tre of the LMC disk and covers a large part of its area, and three 
well-separated smaller-size components. It is reasonable to as¬ 
sume that the extended sources correspond to populations of 
CRs interacting with the ISM, and we explored this possibility 
with a dedicated modelling. 

Under this hypothesis, the emission model was refined by 
adopting a more physical approach to determine the extended 
emission components. In the ~0.I-I00GeV range, the interstel¬ 
lar radiation is dominated by gas-relat ed processes, hadronic in - 
teractions especially (see for instance [Ackermann et al.ll2012bl) . 
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Table 4. Comparison of the emission models 


Model 

logX 

TS 

Ndof 

Background 

-65219.9 

- 

9 

Analytic 

-59167.7 

12104.4 

45 

Emissivity 

-59131.1 

12177.6 

47 

Ionised gas + PSR J0540-6919 

-59462.1 

11515.6 

15 

Ionised gas + 4 point sources 

-59400.0 

11639.8 

27 

All gas + PSR J0540-6919 

-60062.3 

10315.2 

15 

All gas - 1 - 4 point sources 

-59548.9 

11342.0 

27 


From left to right, the columns list the name of the model, 
log(likelihood) values obtained in the model fitting (not including the 
sum of counts term that is model-independent), corresponding test 
statistic with respect to the background model, and the number of de¬ 
grees of freedom in the model. 


The intensity in a given direction can be expressed as the prod¬ 
uct of an emissivity q(Ey), the y-ray emission rate per hydrogen 
atom per unit energy per solid angle (which depends on CR den¬ 
sity and spectrum and on nuclear interaction physics), and a hy¬ 
drogen gas column density (assuming that all gas is pervaded by 
the same CR flux). Instead of searching for a best-fit combina¬ 
tion of 2D Gaussian-shaped intensity distributions, we therefore 
performed an iterative search for a combination of 2D Gaussian- 
shaped emissivity components, which, when multiplied with the 
gas column density distribution of the LMC, provides the best 
fit to the data. We refer to this second model as the emissivity 
model. The gas column density map used in this work includes 
atomic, mole cular, an d ionized hydrogen, and its preparation is 
described in iPaner IL For the dominant atomic and molecular 
hydrogen species, the map is based on ATCA plus Parkes and 
NANTEN data, respectively, assuming optically thin emission 
for atomic hydrogen, and a CO intensity to H 2 column density 
conversion factor of Xco = 7 x 10^°H2cm“^/(Kkms“'). 

In this context, the large-scale component G1 of the analytic 
model is interpreted as arising from a large-scale population of 
CRs spread almost uniformly across the LMC disk and inter¬ 
acting with the gas. An assumption in building the emissivity 
model is that there is a large-scale CR population in the LMC, 
which is described with a 2D Gaussian emissivity distribution 
initially centred on {a, 6) - (80.0°, -68.5°) and having a fixed 
width cr = 3.0° (approximately the angular radius of the galaxy; 
using a flat emissivity profile over the entire extent of the LMC 
was found to degrade the fit by 38 in log -C). We then searched 
for additional emissivity components on top of that. The method 
is very similar to the technique used to build the analytic model: 
at each iteration, a scan over position and size is performed to 
identify a new emissivity component that provides the best fit to 
the data and a significant improvement over the previous itera¬ 
tion. In the fits, a log-parabola spectral shape is assumed for the 
new component because it mimics a local interstellar emission 
spectrum in the ~0.1-lOGeV range fairly well. After all compo¬ 
nents were added, we re-optimised the positions and sizes of the 
components, from the brightest to the faintest in turn, including 
point-like sources (except PI). We also optimised the position of 
the large-scale CR component and found a best-fit position close 
to our initial choice. We then revised the log-parabola spectral 
shape assumption. All components could be described by a sim¬ 
ple power law (see Sect. 15.2b . except for the large-scale compo¬ 
nent, which was found to be more satisfactorily described by a 
log-parabola or a tabulated function derived from the local gas 
emissivity spectrum (see Sect. 15.2b . 



90°00' 84°00' TS^OO' 72°00' 66°00' 

Right Ascension 



Fig. 2. Charts illustrating the layout of the model components. 
Point sources are indicated by green stars. The green circles cor¬ 
respond to the 1-cr extent of the Gaussian emission components 
of the analytic model (top) and to the 1-cr extent of the Gaussian 
emissivity components of the emissivity model (bottom). The 
background map is the total gas column density distribution in 
arbitrary units and a square-root scale. 


As expected, this new approach results in a layout of the 
components which is very close to that of the analytic model. 
There are five extended emissivity components, labelled from EO 
to E4, as listed in Tableland illustrated in Lig.|2 However, two 
of them are adjacent and may well be one single component with 
a shape that is poorly fit by a 2D Gaussian function (components 
El and E3 in the bottom chart of Lig.|3). When we grouped the 
corresponding spatial templates into a single component associ¬ 
ated with a power-law spectral model, we found that the fit was 
hardly degraded (variation in log X < 1). We therefore group the 
components El and E3 of the emissivity model together in the 
following, under the name E1-I-E3. 

One advantage of this method is that the emissivity model 
retains the small-scale structure that the y-ray emission may ac- 
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Table 5. Best-fit spectral parameters of the point-like sources 


Source 

Parameters 

F 200 

0 

0 

txT 

PI (PSRJ0540-6919) 

r= 1.9 ±0.1 
=4.0 ± 1.2 GeV 

(1.0 ± 0.2) X 10-5 

(1.3 ± 0.3) X 10-5 

P2 (PSR J0537-6910 / N 157B) 

r = 2.2 ±0.1 

(0.9 ± 0.2) X 10-5 

(1.1 ±0.2) X 10-5 

P3 (Unassociated) 

r = 2.8 ±0.1 

(5.2 ± 0.6) X 10-5 

(8.9 ± 1.3) X 10-5 

P4 (N 132D) 

r= 1.4 ±0.3 

(1.9 ± 1.4) X 10-5 

(1.9 ± 1.4) X 10-5 


From left to right, the columns list the source identifier used in this work and in parenthesis the established or possible counterpart (see text); 
spectral parameters, a power-law photon index for all sources plus the cutoff energy for PI; and the energy flux in 0.2-100 GeV and extrapolated 
to the 0.1-100 GeV range in MeV cm“^ s“' units. 


tually have, compared to the Gaussian morphologies of the ana¬ 
lytic model. The emissivity model provides a better fit than the 
analytic model, with a difference of 36 in log X for an additional 
2 degrees of freedom, and therefore is the primary model used 
in this work to represent the emission from the LMC. 


3.4. Template fitting 

In iPaner 3 . the emission model providing the highest likelihood 
was a map of the ionised gas column density across the LMC. 
With such templates from observations at other wavelengths, we 
can compare hypotheses about the spatial distribution of the y- 
ray emissio n, which can provide clues about the origin of the 
emission. In lPaner fl. this revealed a surprising correlation of the 
y-ray emission with ionised gas, which traces the population of 
young and massive stars, and a poor correlation with total gas, 
as could have been expected if CRs were almost uniformly dis¬ 
tributed in the LMC. 

We revisited this result with the enlarged data set now avail¬ 
able, and especially in light of the detection of point sources 
within the LMC. We performed hts of the ionised gas and 
total atomic-i-molecularH-ionised g as colum n density distribu¬ 
tions, using the same maps as in iPaner iL For each map, we 
added as other components of the LMC emission model either 
PSR J0540-6919 alone (because it is the only firmly identified 
source), or all four point sources listed in Table[T] with free spec¬ 
tral parameters. The resulting log X are given in Table H) 

The ionised gas template still accounts better for the data 
than the total gas distribution. The impact of the spin temper¬ 
ature used in the preparation of the atomic gas column density 
map is negligible compared to the difference between the two 
templates. But the ionised gas template features a strong peak in 
the 30 Doradus region that pushes sources PI and P2 to negli¬ 
gible TS values (PI in particular is pushed down to a flux level 
below the pulsed fraction of the signal, which is about 75%; see 
lAckermann et al.l[20L5bl) . The ionised gas template is therefore 
inappropriate to model the PSR J0540-6919 point source and 
the extended emission at the same time, while the total gas tem¬ 
plate provides much poorer fits to the data than the emissivity 
or analytic model (differences of ~400 in log X, for ~20 fewer 
degrees of freedom). 


3.5. Summary of all models 

The various emission models tested are summarised in Table 
H] The model yielding the highest likelihood is the emissiv¬ 


ity modefl The corresponding model maps in the 0.2-100 GeV 
range for the background and LMC components are shown in 
the top and middle panels of Fig. [2 respectively, and a residual 
counts map is plotted in the bottom panel of Fig. [2 The distri¬ 
bution of residuals is satisfactory, with a mean close to zero and 
balanced negative and positive residuals. The strongest residual 
excess towards (a, 5) = (79.0°, -72.8°) corresponds to a TS of 
17. 

The outcome of the analysis is insensitive to the size of the 
ROI; we repeated the analysis for larger fields of up to 16° x 16°, 
and the significance of the various emission components as well 
as the normalisations of the isotropic and Galactic diffuse emis¬ 
sion templates were not significantly altered. Since the latter is a 
global model and may be a source of systematics for the specific 
region of the LMC, we also ran an analysis where its different 
components (the various gas templates and the inverse-Compton 
model) were fitted independently; the fit was not significantly 
improved. Overall, the analysis is robust against the background 
model. 

4. Point-like sources 

This section is dedicated to the four point sources found within 
the LMC that are labelled PI, P2, P3, and P4 in Table [T] 
As detailed below, PI is unambiguously identified as pulsar 
PSR J0540-6919 through its characteristic pulsed y-ray emis¬ 
sion. For the other sources, possible associations with known 
objects are discussed, based essentially on their positions and 
spectra. 

Figure 0 ] shows the confidence regions for the localisation 
of point sources P2, P3, and P4 (PI has its position fixed to 
that of the pulsar in optical), plotted on top of a SHASSA Ha 
smoothed mapQ- For each source, we overplotted the positions of 
known SNRs and pulsars in the field, according to the SIMBAD 
database. Figure |2 shows the best-fit spectral shape for each 
source, together with spectral points or upper limits obtained 
from fits in five individual bins per decade in the 0.2-200 GeV 
rang43- Table |2 summarises the best-fit spectral parameters of 
each source and provides the energy flux in the 0.2-100GeV 
band and its extrapolation to the 0.1-100 GeV band. 

^ The significance of the improvement over the analytic model or the 
best template and point sources combi nation cannot be quant ified easily, 
however, since models are not nested dProtassov et al.ll2002f) . 

SHASSA stands for Southern H-Alpha Sky Survey Atlas; more in¬ 
formation at http://amundsen.swarthmore.edu 

^ Spectral points were computed assuming a power-law with fixed 
photon index 2.1 within each bin; upper limits correspond to 95% con¬ 
fidence level and were computed by searching for the flux yielding a 
variation of log X by 1.35 from the optimum. 
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Fig. 3. Top and middle panels: model maps of the background 
components and LMC components of the emissivity model in 
the 0.2-100 GeV band. Black dots indicate the positions of back¬ 
ground sources in the field but outside the LMC boundaries (see 
Sect. l3.1b . Both maps have 0.1° x0.1° pixels and were smoothed 
with a Gaussian kernel with cr = 0.2°. Colours are displayed 
on a square-root scale. Bottom panel: residual counts map in the 
0.2-100 GeV range, after subtracting the fitted emissivity model. 
Colours are displayed on a linear scale. White lines are defined 
in Fig. [1] 


4.1. PI: Pulsar PSR J0540-6919 

Source PI is detected with a TS of 151 and found at a po¬ 
sition that is consistent with that of pulsar PSR J0540-6919 
(also known as PSR B0540-69). PSR J0540-6919 is a 50 ms 
pulsar associated with the ~110 0-year old SNR 0540-69.3 
and PWN N 158A (ISeward et al.l 11984 IWilliams et all 120081: 
iBrantseg et al.l l2014l) . With a spin-down luminosity of £ = 

1.5 X 10^^ erg s ', it has the third largest spin-down power of 
any pulsar currently known, right behind that of the Crab pulsar, 

4.5 X 10^* erg s“' (see the ATNF catalogu^B)- 

The possible contribution of this exceptional pulsar to the 
observed y-r ay signal from the 30 Doradus region was already 
considered in lPaner iL but could not be secured: characteric pul¬ 
sations of the signal were not detected at a significant level, and 
there was some confusion in the morphology of the region re¬ 
garding the presence of a point source in the more extended 
emission of the 30 Doradus region. PSR J0540-6919 is now 
unambiguously detected in GeV y-rays from its characteristic 
pulsations. It is the first y-ray pulsar detected in a galaxy other 
than t he Milky Way and the m ost distant known. A dedicated 
paper (lAckermann et al.ll2015bh reports more details about this 
new y-ray pulsar, such as phase-aligned lightcuryes in different 
bands and a complete spectral energy distribution. 

For consistency with the other sources discussed below, we 
show in Fig.|5]the spectrum of the source obtained from P7REP 
data. The best-fit spectral shape is a power law with exponential 
cutoff, with a statistical significance at the 4cr leyel for the cutoff. 


4.2. P2: Coincident with PSR J0537-6910/N 157B 


Source P2 is detected with a TS of 96 and found at a posi¬ 
tion which is consistent with that of pulsar PSR J0537-6910 
(see Fig. |4|. PSR J0537-6910 is a 16 ms pulsar associated 
with the ~5000-year old SNR 0537.8-6910 and PWN N 157^1 
The whole system presumably results from the explosion of a 
~25 Mq 08-09 star and is likely associated with the O B as¬ 
sociation LH99 (IChen et alJ 12004 iMicelotta et akl l2009l). The 
pulsa r was discoyered from X-ray obseryations ( Marshall et al.l 
Il998h and is the most rapidly spinning and the most power¬ 
ful young pulsar known, with a spin-down luminosity of £ = 
4.9 X 10^® erg s''. 

PWN N 157B was detecte d in TeV y-rays with H.E.S.S. 
(lAbramowski et a02012LlT015h . It is the first and only PWN de¬ 
tected outside of the Milky Way at these photon energies, which 
can be explained by the enormous spin-down power of the pulsar 
combined with a rich infrared photon field for inyerse-Compton 
scattering. The emission in the £emu-LAT range is obseryed at 
a leyel of ~ 2-3 xlO'*^ ergcm'^ s'*, quite similar to that of the 
TeV emission detected with H.E.S.S.. 

As shown in Fig. |5] howeyer, the spectrum of the source 
is rather flat oyer the 0.5-100GeV range, with a photon in¬ 
dex of 2.2+0.1. The origin of the emission is unclear. It may 
come from the pulsar, the SNR, the PWN, or any combination 
of those in unknown proportions. The nearly flat or eyen de¬ 
clining spectrum seems hard to reconcile with emission from 
a PWN alone: in a typical scenario, emission from a PWN in 
the Fermi-LAiT range is expected to haye a hard spectrum as it 


® See lManchester et ^ l l2005h and the ATNF Pulsar Catalog version 
1.52 at http://www.atnf.csiro.au/people/pulsar/psrcat/ 

^ The complete nomenclature is LHA 120-N 157B hut we use 
N 157B for short. The same remark applies to N 158A, N 132D, Nil, 
and other sources discussed further down. 
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probes the low-energy side of an inverse-Compton hump peak¬ 
ing at ~T eV energies. In tha t respe ct, the spectral model con¬ 
sidered in lAbramowski et all (1201 5h for the TeV data does not 
match our GeV observations at all. This situation is somewhat 
reminiscent of the Vela X PWN, for which GeV and TeV emis¬ 
sions at similar levels were postulated to arise from two d iffer- 
ent populations of emitting particles (lAbdo et al.ll2010l3) . The 
observed emission may also result from the superposition of 
a pulsar’s curved spectrum peaking at ~lGeV and a PWN’s 
spectrum ramping up and dominating at ~ 10 GeV The detec¬ 
tion of pulsations from the source and an analysis of its off-pulse 
emission would have been key to establish such a scenario. We 
searched for pulsa tions in the y-ray signal in the same way as for 
PSR J0540-6919 (lAckermann et al.l2015bl) . using an ephemeris 
for PSR J0537-6910 based on RXTE observations covering the 
first 3.5 years of the Fermi mission. No evidence for y-ray pul¬ 
sations was found, with a pulsation significance below Icr. On 
longer time scales, the source shows no evidence of variability 
down to a monthly basis. 

With a total luminosity of 5.2 x lO^^ergs^' in the 0.1- 
100 GeV range, this source is remarkable. Further studies in¬ 
volving the complete picture of the plerionic system N 157B 
are required to assess the possible origin of the GeV emission. It 
is clear, however, that the accuracy of the localisation allows us 
to exclude superbubble 30 Doradus C as the main contributor to 
the emission (s ee Fig.lHl. The latter was detected in TeV y-rays 
with H.F.S.S. (lAbramowski et alll2015h and might have had a 
GeV counterpart detectable with Termi-LAT. An upper limit on 
the emission from this object is given below, in Sect. 14.61 

4.3. P3: Unassociated source 

Source P3 is detected with aTS of 116 and found at a position 
which is consistent with that of HII regions NGC 2029/NGC 
2032. Supernova remnant DFM L 241 and Seyfert 1 galaxy 2F 
1445 are also in the field but they lie at a greater distance from 
the best-fit position, at the edge of the 2cr confidence region. 

The source shows no evidence of variability down to a 
monthly basis. The most notable feature of this source is its very 
soft power-law spectrum, with a photon index of 2.8 + 0.1. If 
source P3 is a background active galactic nucleus, such a soft 
spectrum tends to f avour a flat spectrum rad io quasar type over 
a BL Lacertae type (lAckermann et al.ll201 lah . 

4.4. P4: Coincident with SNR N 132D 

Source P4 is detected with aTS of 28 and found at a position 
which is marginally consistent with that of SNR N 132D, just 
outside the Icr confidence region (see Fig. Ell. Supernova rem¬ 
nant N 132D belongs to the rare class of O-rich remnants, to¬ 
gether with objects such as Cassiopeia A, Pup pis A, or SNR 
0540 -69.3, and it has a kinematic age of 2500 yr Jvogt & Donital 
|201 iL and references therein). Its chemical composition indi¬ 
cates that N 132D is the remnant of the core-collapse supernova 
of a progenitor sta r with initial mass in the ~35-75Mo range 
(iFrance et alj|2009l) . The shock wave in the remnant is interact¬ 
ing with surrounding mater ial, and N 132D is the b rightest rem¬ 
nant of the LMC in X-rays (iBorkowski et al.ll2007h . 

The association of our weak GeV source with N 132D is 
supported by the recent marginal detection using H.F.S.S. of a 
TeV source spatially coincident with N 132D (with a sub-arcmin 
localisation accuracy compare d to the few arcmin allo wed by 
Fermi-LAST for this source; see lAbramowski et al.ll2015]) . In the 


Fenni-LAiT range, the source has a hard spectrum with photon 
index 1.4 ± 0.3 and a differential flux of ~ lO^^MeV cm“^ s“' 
at 100GeV. In the H.F.S.S. band, the differential flux is ~ 2 x 
10“^ MeV cm“^ s“' at 1 TeV, with a photon index 2.4+0.3. If the 
y-ray source is N 132D, its spectral energy distribution in the y- 
ray domain may have a peak around 100 GeV, possibly a broad 
one given the uncertainty on the spectral index in the 1-100 GeV 
range. The spectrum presente d here is a factor ~ 5 ab ove the 
hadronic model considered in lAbramowski et akl (1201 5l) . which 
may be accommodated by a modest tuning of the model param¬ 
eters; in contrast, the leptonic model seems excluded because of 
a much lower predicted flux and possibly also an inconsistent 
spectral slope. 

4.5. Systematic uncertainties 

The properties of the four point sources discussed above may 
be affected by systematic errors, introduced by selecting a par¬ 
ticular method to build the complete LMC emission model. As 
exposed in Sect. |2 two models provide a good fit to the data: 
the emissivity model, which yields the highest likelihood and is 
the reference one, and the analytic model, which gives a slightly 
lower likelihood but can be considered as a competitive alterna¬ 
tive. The other models tested in this work were far from equaling 
these two (see TableEll- 

Interestingly, the 2D Gaussian extended components of the 
analytic model have no fine structures compared to the extended 
components of the emissivity model, which were built from the 
gas column density distribution. The analytic model can there¬ 
fore be used to test whether such fine structures can bias the 
inferred properties of point sources. 

When using the analytic model, the spectral parameters of 
the four point sources within the LMC turned out to be consistent 
within Icr statistical uncertainties with those obtained with the 
emissivity model. The differences in the mean values of photon 
and energy fluxes of point sources between the two models were 
found to be < 50% of the statistical uncertainties. 

The inferred properties of the sources may also be affected 
by uncertainti es in the characte r isation of the LAT instrument. 
According to lAckermann et al.l (1201 2al) . the magnitude of the 
systematic error arising from uncertainties in the effective area, 
point spread function, and energy dispersion, is about 0.1 on 
the spectral index and about 10% on the photon and energy 
fluxes (for the analysis of a not-too-hard point source in the 0.1- 
100 GeV range). These systematic uncertainties are therefore of 
the same order as the statistical uncertainties on spectral indices 
and fluxes for the three brightest point sources in the LMC. 

4.6. Upper limits on selected objects 

Some particularly interesting LMC objects that were or could 
have been expected to be y-ray sources were undetected in the 
Fermi-LAT observations we analysed. This is the case of su¬ 
perbubble 30 Doradus C (a rare superbubble with non-thermal 
emission, recently detec t ed at TeV energies wit h H.F.S.S.; see 
lAbramowski et al.ll2015l: IWmaguchi et al.ll201(jh . or SN 1987A 
(the closest supernova in modern times, predicted to be accel- 
erating particles and to be increasingly shining in y-rays; see 
iBerezhko et al.l 1^1 Ih . Another possible source of GeV emis¬ 
sion was N 11, the second most active star-forming region of 
the LMC after 30 Doradus, and maybe an o lder version of it 
(iRosado et alJ[l9^ IWalborn & Parked [l9^ . The LMC also 
harbours many SNRs that are very well studied in X-rays and at 


8 









































LAT collaboration: A deep view of the Large Magellanic Cloud 


radio wavelengths, and for which constraints on the y-ray emis¬ 
sion might be useful. 



Right Ascension 



82 ° 00 ' 40 ' 20 ' 81 ° 00 ' 80 ° 40 ' 

Right Ascension 


Fig. 4. Localization maps for point sources P2, P3, and P4 (from 
top to bottom; PI is identified as PSR J0540-6919 and thus not 
shown). Green contours correspond to 1,2, and 3cr confidence 
regions derived using the emissivity model (4 and 5cr confidence 
regions were added to P2 to emphasise that it is not consistent 
with 30 Dor C); the cyan dots correspond to the best-fit positions 
derived using the analytic model. The background is a SHASSA 
Ha smoothed map on a log scale. Overlaid as red dots and cir¬ 
cles are the position of known SNRs, including the plerions as¬ 
sociated with PSRs J0540-6919 and J0537-6910, together with 
other noteworthy objects such as SN1987A or 30 Doradus C; 
magenta pluses indicate the positions of known pulsars, whil| 
yellow crosses indicate the positions of known active galactic 
nuclei. 
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Table 6. Upper limits on undetected sources 


Source 

^ = 2.0 

.? = 2.5 

SN 1987A 

7.8 X 10^’ 

6.6 X 10-’ 

B0505-67.9 (DEML71) 

5.9 X 10-^ 

6.1 X 10“’ 

B0506-68.0 (N 23) 

6.7 X 10^^ 

7.0 X 10“’ 

B0509-68.7 

3.2 X 10-^ 

3.8 X 10-’ 

B0509-67.5 

1.3 X 10-^ 

1.4 X 10-’ 

B0519-69.0 

5.6 X 10-^ 

5.9 X 10-’ 

B0525-66.0 (N 49B) 

3.1 X 10^^ 

3.9 X 10-’ 

B0525-66.1 (N49) 

3.1 X 10^^ 

3.9 X 10-’ 

B0534-69.9 

3.2 X 10“^ 

2.7 X 10“’ 

B0535-66.0 (N 63A) 

4.8 X 10-^ 

4.8 X 10-’ 

B0548-70.4 

1.6 X 10-^ 

2.0 X 10-’ 

30 Doradus C 

1.1 X 10 

8.9 X 10“’ 

N 11 

7.1 X 10^^ 

7.1 X 10“’ 


Upper limits correspond to a 95% confidence level and were com¬ 
puted under the assumption of a power-law spectrum. Column 1 lists 
the sources for which upper limits were computed, and Cols. 2 and 3 
give the flux limits in units of MeV cm“^ s“' over the 1-10 GeV interval 
for two power-law photon indices. All SNRs were modelled as point 
sources, while Nil and 30 Dor C were modelled as 2D Gaussian in¬ 
tensity distributions with cr = 0.4° and cr = 0.1° , respectively. 

Table 7. Comparison of spectral models for component EO 


Model 

logX 

Ndof 

Power law 

-102523.4 

2 

Power law with cutoff 

-102499.6 

3 

Broken power law 

-102497.3 

4 

Log parabola 

-102496.3 

3 

Tabulated 

-102495.7 

1 


logX values correspond to a fit over 0.2-10 GeV with 17 logarith¬ 
mic bins (because component EO is not significantly detected above 
10 GeV), hence the difference with values in Table|4] 

Based on the emissivity model described previously, we de¬ 
rived upper limits on the 1-10 GeV flux from the three objects 
mentioned above and the ten brightest remnants in X-rays, ac¬ 
cording to the MCSNR databas^!]. These are summarised in 
Table |6] under two assumptions for the photon index over 1- 
10 GeV. 

For SN1987A, our upper limit on the 1-10 GeV flux is 
not constraining whe n compared to the predictions made by 
iBerezhko et al.l (1201 ih : our limit is a factor of 3 above their pre¬ 
dicted flux for 2010 and remains above the maximum flux pre¬ 
dicted for 2030. For 30 Dor C, our upper limit is not constraining 
ei ther when compared to th e y-ray emission models considered 
in lAbramowski et^ (l2015h : our limit is a factor of 6 above the 
most optimistic prediction of the hadronic model. 


www.mcsnr.org 






Fig. 5. Spectra of the four point-like sources found in the LMC, 
from top to bottom PI to P4, as listed in Table[T] We show spec¬ 
tral points obtained by fits in individual bins, together with the 
best-fit spectral model from a binned likelihood analysis over 
the full energy range. The spectrum of the weakest source P4 
features coarser binning in energy. Upper limits (arrows) corre¬ 
spond to a 95% confidence level. 
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5. Extended sources 

This section is dedicated to the extended sources found in the 
modelling of the y-ray emission from the LMC. We used as ref¬ 
erence the emissivity model and the components labelled EO, 
El, E2, E3, and E4 in Tabled We recall that components El 
and E3 can be grouped in the fit of the model, because leaving 
them independent does not provide a significant improvement; 
however, they need to be considered separately in some parts of 
the discussion below. 

We first discuss the localisation and spectrum of each ex¬ 
tended source. Eigure|6]illustrates the layout of each region asso¬ 
ciated with extended y-ray emission. EigureQ shows the best-fit 
spectral shape for each component, together with spectral points 
and upper limits derived over the 0.2-200 GeV range 0. Table 
[8] summarises the best-fit spectral parameters of each extended 
emission component and provides the energy flux in the 0.2- 
100 GeV band and its extrapolation to the 0.1-100 GeV band. 

5.1. Localisation 

The large-scale component labelled EO is by construction emis¬ 
sion coming from the whole LMC disk. A large part of it prob¬ 
ably results from a large-scale population of CRs accumulated 
over long time scales and interacting with interstellar gas (see 
Sect. I5.3b . The origin of the emission from the three to four 
smaller-size extended components is more uncertain. It is not 
correlated with prominent features of the gas column density 
distribution, which suggests at least two possibilities: these rela¬ 
tively gas-poor regions harbour an overdensity of CRs as a result 
of a local concentration of CR sources, or they host a population 
of many discrete y-ray emitters with individual fluxes below the 
detection threshold. To explore both options, we searched for 
correlation with structures or objects in the LMC. In Eig.|6] the 
location and size of emissivity components (from the emissivity 
model) and emission components (from the analytic model) are 
plotted on top of a SHASSA Ha smoothed map, which indirectly 
reveals the distribution of young stars. We also overlaid the dis¬ 
tributions of various objects: pulsars, which are the dominant 
class of GeV sources in the Galaxy; Wolf-Rayet stars, which re¬ 
veal where the most massive and shortest-lived stars are located; 
high-mass X-ray binaries (HMXBs); stars of spectral type from 
BO to B3, which offer a longer perspective on recent star forma¬ 
tion because their typical lifetime is a few lOMyr; supergiant 
shells (approximat ed as circles), w hich are structures with typi¬ 
cal ages 5-lOMyr (iKim et al.lll999h . testifying to strong releases 
of kinetic energy in the ISM, from which CRs could be acceler¬ 
ated and confined. 

No obvious correlation appears to be valid for all three re¬ 
gions at the same time. The most striking feature is the coinci¬ 
dence of emissivity component E2 with a large void delimited 
by a supergiant shell. The cavity is relatively empty of objects, 
at least those considered above; young stars of the Wolf-Rayet 
class to early B spectral type tend to be found on the periph¬ 
ery (the cluster of B stars on the southe rn edge of the cavity is 
NGC 2004; see iNiederhofer et al.ll2015h . This also seems to be 
the case for emissivity component E4, although less remarkably 
so. The situation is less clear for emissivity components El and 
E3 owing to the crowding of the region. Most young stars are 
still found outside the 1-cr contour of the emissivity component; 
in particular, the rich 30 Doradus area and the surrounding HII 
regions are clearly offset from the core of the emissivity compo- 


® Data points and upper limits computed as described in footnote 5. 


nent. Maybe the most interesting feature is the correlation with 
the highest concentration of supergiant shells in the LMC, which 
testifies to intense star formation activity over the past 5-10 Myr. 
The absence of conspicuous voids associated with some of these 
supergiant shells may also indicate that projection effects com¬ 
plicate such a search for spatial correlation. 

5.2. Spectra 

The spectrum of the large-scale disk component labelled EO 
shows some curvature (see Eig. |7]i. We searched for the best 
model for it among a power law with exponential cutoff, a bro¬ 
ken power law, a log-parabola function, and a tabulated func¬ 
tion derived from the local gas emissivity spectrum (as a result 
of its interactions with CRs). Results are summarised in Table 
|7] The best analytical model is the log-parabola. The broken 
power law is almost as good (logX decreases by 1 despite 1 
additional degree of freedom), while the power law with cut¬ 
off is disfavoured (log -C decreases by 3 for the same number 
of degrees of freedom). The significance of the curvature al¬ 
lowed by the log-parabola model compared to the power law is 
7.4cr. The tabulated function model slightly improves upon the 
log-parabola, with a likelihood increase by about 1 for 2 fewer 
degrees of freedom. At this stage, however, we avoided being 
too specific regarding the emission associated with component 
EO. We therefore used the best-fit log-parabola function, with 
spectral parameters listed in Table [8] and corresponding spec¬ 
trum shown in Eig. |7] An interpretation in terms of CRs inter¬ 
acting with the gas is presented below, among other possibili¬ 
ties, and the use of the tabulated function as a spectral model 
is discussed there (the similarity between the log-parabola and 
tabulated function models can be inferred from Eig. Ell. 

The spectra of components E1-I-E3, E2, and E4 can be de¬ 
scribed with simple power laws having relatively hard photon 
indices of ~2.0-2.2 (see Eig. |7]i. We evaluated whether these 
small-scale components can be distinguished from large-scale 
component EO in terms of spectral properties. At low energies 
< 2 GeV, they have spectra consistent with that of EO (a power 
law with photon index 1.9, as inferred from the broken power- 
law fit). At energies > 10 GeV, only E2 and E4 are significantly 
detected (with TS > 20), while E1-I-E3 and EO are no{3- In a 
binned likelihood fit over the entire energy range, replacing the 
power-law assumption for components E1-I-E3, E2, and E4 by 
the tabulated function model tested on component EO yields the 
following results: the fit is hardly degraded for E4, with vari¬ 
ations in log X. that are smaller than the number of degrees of 
freedom lost; in contrast, the fit for components E1-I-E3 and E2 
degrades by about 10 and 20 in log X, respectively, for 1 degree 
of freedom lost. These results might seem to contradict the fact 
that E4 is significantly detected above 10 GeV, while E1-I-E3 is 
not. The spectra of these two components appear somewhat ir¬ 
regular in Eig. |7] when compared for instance to those of point 
sources (for E4, this is more obvious when using a finer energy 
binning). A possible explanation is that these extended compo¬ 
nents encompass different sources with different spectra (as con¬ 
sidered in Sect. 15.4b . such that a single power law or the tabulated 
model are just first-order approximations. 


That components E2 and E4 could be detected above 10 GeV at lev¬ 
els ~ 3 x 10“® MeV cm“^ s“‘ while component EO is only constrained 
by an upper limit at a similar magnitude is due to their smaller size; this 
results in higher surface brightnesses, such that the smaller-scale com¬ 
ponents can be detected against the diffuse background at comparatively 
lower fluxes. 


11 








LAT collaboration: A deep view of the Large Magellanic Cloud 


Table 8. Best-fit spectral parameters of the extended components 


Source 

Parameters 

E200 

Fioo 

EO (large-scale disk) 

a = 2.17 ±0.04 
b = 0.22 ± 0.02 
Et = 1.0 GeV 

(5.9 ± 0.3) X 10-5 

(6.7 ± 0.4) X 10-5 

E1+E3 (west of 30 Doradus) 

r = 2.13 ±0.05 

(2.1 ±0.2) X 10-5 

(2.5 ± 0.2) X 10-5 

E2 (northern LMC) 

r = 2.00 ± 0.06 

(1.3 ± 0.2) X 10-5 

(1.4 ± 0.2) X 10-5 

E4 (western LMC) 

r = 2.13 ±0.07 

(0.9 ± 0.1) X 10-5 

(1.1 ±0.1) X 10-5 


Note to the Table: From left to right, columns are: the source identifier used in this work; the spectral parameters, a power-law photon index for the 
last three sources, and log-parabola parameters for EO; the energy flux in 0.2-100 GeV and extrapolated to the 0.1-100 GeV range in MeV cm“^ s“* 
units. 


The total energy flux from the extended emission compo¬ 
nents in our model amounts to 1.2 x 10“"* MeVcm^^ s“' over 
the O.l-lOOGeV range. In Fig. 0 the total spectrum of that 
emission is plotted against the spectra of the point sources 
for which we have an established (PSR J0540-6919) or likely 
(PSR J0537-6910 and N 132D) association with an LMC ob¬ 
ject. To build this total spectrum, we added the best-fit spec¬ 
tral functions of all extended components together. Emission of 
extended nature clearly dominates the GeV emission from the 
LMC. It is interesting to note, however, that the summed emis¬ 
sion from only two discrete sources (PSR J0540-6919 and the 
source coincident with PSR J0537-6910) amounts to about one 
third of the extended emission flux at 100 MeV, and to about 
20% of it at 10 GeV, and that the total point-like emission has a 
similar spectrum to the total extended emission. 

5.3. Cosmic-ray population 

The emissivity model introduced above to describe the extended 
emission from the LMC provides direct estimates of the emissiv¬ 
ity distribution over the galaxy. Before addressing this quantita¬ 
tively, however, we wish to emphasise that this is an assumption 
about the nature of the extended emission we observe: it predom¬ 
inantly results from hadronic interactions between CR nuclei and 
interstellar gas. The increase in likelihood it leads to, compared 
to the more phenomenological description of the analytic model, 
seemingly supports that interpretation, but this is not a formal 
proof that the latter is valid. There are other possible contribu¬ 
tions to the extended emission, such as unresolved populations 
of discrete sources such as pulsars, or inverse-Compton emission 
from CR electrons. They were not tested against the data because 
their modelling is much more uncertain, but they are addressed 
below and in Sect. 15.41 Lor this reason, the emissivity estimates 
given below should be considered as upper limits. 

Lor each extended emission component of the emissivity 
model, the gas emissivity spectrum is obtained by dividing the 
fitted spectrum for that component by the average column den¬ 
sity over the region, weighted by the 2D Gaussian emissivity 
profile. The resulting emissivity spectra are plotted in Lig.|9] Lor 
comparison, we plotted the loca l Galactic emissivity spectrum 
determined bv ICasandiianI (1201 2h from high-latitude Fermi-hAJ' 
observations in the ~50MeV-40GeV range. 

In terms of spectral shape, the emissivity spectra are not all 
consistent with that measured locally in the Milky Way. As men¬ 
tioned in Sect. 15.21 component E4 has a spectrum consistent with 
the local emissivity, while E1-I-E3 and E2 significantly differ 
from it. Splitting E1-I-E3 into independent components El and 
E3 reveals that El is consistent with the local emissivity spec¬ 


trum, while E3 is not. That the emissivity spectral shapes of El 
and E4 are consistent with the local one suggests that the spec¬ 
tral distribution of CRs is similar in these regions and in the local 
medium, at least up to CR energies of a few hundred GeV (re¬ 
sponsible for emission at a few tens of GeV). In these conditions, 
the differences in the normalisation of the emissivity can be re¬ 
lated directly to a difference in CR density, as discussed below. 

In terms of normalisation, the emissivity components depart 
from the local value. We emphasise that emissivity components 
were modelled with 2D Gaussian spatial profiles, which implies 
a decrease of the emissivity away from the best-fit central posi¬ 
tion; in the following, the emissivities of the various components 
are given as peak values at the centre of the Gaussian profiles. 
The large-scale disk peak emissivity is 33% of the local emissiv¬ 
ity. As mentioned in Sect. l3.3l this component can be interpreted 
as resulting from CRs accumulated in the LMC over time scales 
that are long compared to the recurrence time of CR injection 
in the ISM (by supernovae). The difference in emissivity level 
means that the large-scale density of CRs is about one third of 
the local Galactic density in the inner regions of the LMC and 
decreases to about 10% in the outskirts (by construction, because 
the size of the component was fixed). At 1 GeV, the emissivities 
of components El and E4 are 40% and 60% of the local emissiv¬ 
ity, respectively. Although this implies an underdensity in CRs 
with respect to the local medium, it means an enhancement in 
CR density with respect to the large-scale peak level in the LMC 
by a factor 3 for component E4. Components E2 and E3 have 
an emissivity that is twice the local value at 1 GeV, or about 6 
times the large-scale peak level in the LMC. Because of the sig¬ 
nificantly harder emissivity spectrum of the E2 component, the 
enhancement in CR density increases with energy, up to about 8 
times the local value at 10 GeV or 24 times the peak large-scale 
level in the LMC. The main source of uncertainty on these emis¬ 
sivity estimates is the gas content of the LMC (see Appendix 
A). 
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Fig. 6. Charts illustrating the layout of each region associated 
with small-scale extended y-ray emission. From top to bottom: 
regions around emissivity components E1-I-E3, E2, and E4. The 
background is a SHASSA Ha smoothed map on a log scale. The 
green circles are the 1-cr extent of the Gaussian emissivity com¬ 
ponents of the emissivity model, the cyan circles are the 1-cr ex¬ 
tent of the Gaussian emission components of the analytic model. 
Overlaid are pulsars (magenta pluses), Wolf-Rayet stars (white 
diamonds), stars of spectral type BO-3 (blue dots), HMXBs (red 
dots), and supergiant shells (yellow circles). 









H 

J 





El+EB (west O' 

30 Doradus) 
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Fig. 7. Same as Fig. |5] for the four extended emission compo¬ 
nents, from top to bottom EO, E1-I-E3, E2, and E4, as listed in 
Table [3] The spectrum of EO has five bins per decade and a sin¬ 
gle upper limit above 20 GeV computed for a fixed power-law 
photon index of 2.7; the spectra of other components have four 
bins per decade. 
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Fig. 8. Comparison of the total emission of extended nature with 
that of point sources for which an established or likely associa¬ 
tion with an LMC object exists (as indicated in the legend). Best- 
fit analytic spectra were used and summed over components EO 
to E4 to produce the black curve for the diffuse contribution. 


Fig. 9. Emissivity spectra of the five components of the emissiv- 
ity model, obtained by dividing the fitted spectrum for that com¬ 
ponent by the average column density over the region, weighted 
by the 2D Gaussian emissivity profile. The dashed line is the lo¬ 
cal Galactic emissivity spectrum from ICasandiianI (l2012h . given 
here for comparison. 


Interestingly, components whose spectra are compatible with 
the local emissivity are those with lower normalisations, while 
components that are significantly harder are associated with 
higher normalisations. The high CR densities and harder emis¬ 
sivity spectra in components E2 and E3 suggest an accumulation 
of fresh CRs in these regions; these CRs would have been re¬ 
leased relatively recently by a spatial and temporal concentration 
of sources, and they would not have had the time to diffuse away 
substantially, which preserves their hard injection spectrum (re¬ 
flected in the emissivity spectrum). This situation, and the possi¬ 
ble association of emission components with cavities and super¬ 
giant shells (see Sect. EB, is re miniscent of the Cygnus “co¬ 
coon” found in the Milky Way (lAckermann et al.ll201 Ibl) . al¬ 
though on much larger spatial scales. But it is difficult to firmly 
establish such a scenario, especially given the lack of clear spa¬ 
tial correlation with other objects (see Sect. EB. Overall, the 
CR population inferred above corresponds to 1.4 x 10^^ erg in 
the form of ~l-100GeV CRs in the disk, with an additional 
0.9 X 10^^ erg contained in the smaller-scale regions; given the 
supernova rate in the LMC, this implies a ~ 1 Myr residence time 
for CRs in the LMC disk (see Appendix B). 

An alternative explanation for the relatively hard spectra of 
extended components is that inverse-Compton emission domi¬ 
nates the signal from these sources. In models of the steady- 
stat e CR-induced interstellar emission from star-forming galax¬ 
ies (lMartinll2014l) . the inverse-Compton spectrum is harder than 
the pion decay spectrum above a few lOOMeV and up to a few 
lOOGeV, and it can account for a significant fraction ~20% 
or mo re of the O.l-lOOGeV flux. Erom the models of iMartiiil 
(l2014l) . a small synthetic galaxy was found to be a possible 
representation of the LMCFH. It yields a total interstellar dif¬ 
ferential flux of 2.5 X 10“^ MeV cm“^ s * at IGeV, consistent 
with the observed value (see Eig.[8]l. The corresponding inverse- 
Compton emission spectrum was found to fit the Termi-LAT 
spectra of components E2 and E4 as well as power laws, and 
even to slightly improve the fit for component E1-I-E3 (by about 
2 in log for 1 fewer degree of freedom); in contrast, the spec¬ 

* * The model is an atomic gas disk of 2 kpc radius with a molecular 
core of 100 pc and a density n = 20 H2 cm“^. 


trum of the large-scale EO component is poorly represented by 
the inverse-Compton spectrum, with a degradation of 62 in in 
log X., for 2 fewer degrees of freedom. Erom the spectral point 
of view, the small-scale extended components are therefore con¬ 
sistent with an inverse-Compton origin of the emission. Erom 
the spatial point of view, however, the locations of these com¬ 
ponents are somewhat unexpected, at least for components E2 
and E4. These are found in regions of the LMC that are rela¬ 
tively far away from strong sources of radiation, infrared or op¬ 
tical, and that are not known for particularly high strengths of 
the interstellar radiation field ([Bernard et alj|200^ . In the ab¬ 
sence of particularly abundant targets for the inverse-Compton 
scattering, emission from these regions would therefore require 
localised enhancements of the CR electron density, similar to 
the CR nuclei enhancement considered previously. We then face 
the same problem that these regions are not clearly correlated 
with particularly high concentrations of massive stars or SNRs 
or any other possible tracer of the past CR injection activity (see 
Eig.ig. Independent of this consideration, such accumulations of 
high-energy electrons may be expected to be clearly detectable 
in radio synchrotron emission (provided the magnetic field is not 
anom alously low), but nothing particular is observed (iMao et al.l 
l2012h . The case of component E1-I-E3 seems different. Being 
closer to the 30 Doradus region and also near the so-called opti¬ 
cal bar, the situation is more favourable regarding the availability 
of high-energy electrons and target photons. 

5.4. Unresolved population of sources 

As suggested earlier, extended emission from the LMC may 
arise from an unresolved population of low-luminosity y-ray 
sources instead of CRs interacting with the ISM. This may be 
especially true for emission from the small-sized regions (G2 to 
G4 in the analytic model, El to E4 in the emissivity model). The 
number of such objects does not need to be large. As an exam¬ 
ple, model component E2 for extended emission in the northern 
LMC should be replaced by at least five point-like sources with 
TS ~25-50 to provide a similar fit likelihood (but the number of 
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degrees of freedom involved in searching for the best-fit prop¬ 
erties of these point-like sources is four times that needed for a 
single extended component). Components El to E4 have simi¬ 
lar fluxes, so >20 objects would be needed to account for their 
cumulated emission. 

Since the dominant class of identified GeV sources in the 
Milky Way is pulsars (lAcero et al.l 1201 5h . it is reasonable to 
assume that two dozen y-ray pulsars with luminosities much 
lower than that of the exceptional PSR J0540-6919 could ap¬ 
pear as extended emission from specific regions of the EMC. 
Interestingly, localisation maps in Eig. 0 already show a hand¬ 
ful of known pulsars that coincide with the extended emission 
components: 11 out of 25 pulsars are found within the 1-cr an¬ 
gular radius of regions G2 to G4. Many others may lie there and 
be currently undetected in radio and X-ray surveys (an estimated 
number of potentially observable pulsars in the EMC is of the or¬ 
der of several 10^; see lRidlev et al.ll20T^ . Some of the extended 
emission components have hard spectra (for instance E2), how¬ 
ever, with significant emission above 10 GeV; this seems difficult 
to reconcile with the now established fact that the majority of 
middle-aged pulsa rs have spectra cutting off beyond a few GeV 
(lAbdo et al.ll2013h . On the other hand, pulsars tend to be associ¬ 
ated with PWN, which are expected to have hard spectra at GeV 
energies, and two out of the four sources found in the EMC are 
detected beyond 10 GeV. 

An unresolved population of GeV sources may also include 
low surface-brightness SNRs. The y-ray luminosities and spec¬ 
tra of SNRs likely depend on many factors, including type of the 
progenitor, age of the remnant, and environmental conditions. In 
the Galaxy, the integrated emission from SNRs was put forward 
as a possible explanation for the observed harder 1-100 GeV 
emission from the inner regions, compared to expectations based 
on la rge-scale interstellar radiation alone (I Volk & Berezhkol 
l2013h . In the EMC, SNRs could therefore add their contribu¬ 
tion to that of PWNe above 10 GeV, where pulsar emission is 
expected to rapidly decrease, and thus account for the hard emis¬ 
sion from some extended emission components. But SNRs and 
PWNe are relatively short-lived, with typical lifetimes of order 
10"^ yr, and this is difficult to reconcile with the observation that 
some extended emission components seem to coincide with su¬ 
pergiant shells with ages of several Myr. 

Overall, both unresolved populations of sources and CR- 
induced interstellar emission are possible explanations for at 
least some fraction of the extended y-ray emission from the 
EMC. The two options are difficult to distinguish based on GeV 
observations alone, and considering observations at other ener¬ 
gies will most likely be the key to clarifying the situation. Eor in¬ 
stance, future ~TeV observations with the Cherenkov Telescope 
Array may determine whether the fraction of the emission is 
indeed diffuse above 100 GeV by detecting extended emission 
and/or by finding new discrete sources. At the other end of the 
electromagnetic spectrum, Planck will soon provide a new view 
of the synchrotron emission from the EMC, thus revealing how 
the leptonic component of CRs is distributed, and future ~GHz 
observations with SKA may refine the picture even more and 
increase the number of known pulsars in the EMC. 

6. Conclusion 

We have analysed Fermi-LP^ observations of the EMC over 
about 73 months of operation. The analysis was made over 
the 0.2-100 GeV range. After subtracting foreground and back¬ 
ground signals, the y-ray emission from within the boundaries 
of the EMC can be modelled as a combination of four point-like 


objects, a large-scale component about the size of the galaxy, and 
three to four small-scale regions about 1-2° in angular extent. 

One of the point sources is unambiguously identified as pul¬ 
sar PSR J0540-6919 through its characteristic pulsations; it is 
the first y-ray pulsar detected outside the Milky Way and has 
the highest isotropic y-ray luminosity. Another source is spa¬ 
tially coincident with plerion N 157B and its powerful pulsar 
PSR J0537-6910; no pulsations from the latter are detected, 
and the relatively hard spectrum extending up to >50 GeV sug¬ 
gests that the emission is not contributed by the pulsar alone, but 
may also be related to the TeV source detected with H.E.S.S.. 
A third point source is spatially coincident with the bright SNR 
N 132D, also detected at TeV energies by H.E.S.S.; this weak 
source is detected above 10 GeV with a relatively hard spectrum 
that implies a broad peak in the spectral energy distribution at 
~ 100 GeV. A fourth point source is currently unidentified and 
could be a background active galactic nucleus; it shows no sign 
of variability and its spectrum is soft with a photon index of 
~ 2 . 8 . 

The prominent 30 Doradus massive star-forming region 
is quite bright in GeV y-rays mainly because of the pres¬ 
ence of PSR J0540-6919 and the source coincident with 
PSR J0537-6910. In contrast, the second most active star¬ 
forming region in the EMC, N 11, does not show a remark¬ 
able enhancement of y-ray emission on top of the large-scale ex¬ 
tended emission. SN1987A and 30 Doradus C are not detected. 

In spite of resolving individual GeV sources in the EMC for 
the first time, emission of extended nature is still significantly 
detected and accounts for the majority of the EMC y-ray flux. 
The nature of this extended emission remains unclear. In par¬ 
ticular, the small-scale extended regions do not correlate with 
known objects or structures in the EMC, except perhaps having 
some relation to cavities and supergiant shells. 

A model in which the extended emission is interpreted as in¬ 
teractions of a spatially inhomogeneous population of CRs with 
interstellar gas provides a satisfactory fit to the data. Under this 
assumption, the EMC would be filled with a large-scale popu¬ 
lation of CRs with a central peak density of one third the local 
Galactic density; in addition, the small-scale extended emission 
components would be regions where the density of CRs is en¬ 
hanced by factors of 2 to 6 at least. The y-ray emission from 
some of these regions features relatively hard spectra, which can 
be interpreted as resulting from a more energetic and possibly 
younger population of CRs compared to that filling the EMC on 
large scales. An alternative scenario for the origin of the small- 
scale extended emission is that of an unresolved population of 
sources. Two dozen or more weak point sources can mimic the 
extended nature of the emission. Currently unknown pulsars and 
their associated nebulae or supernova remnants would be obvi¬ 
ous candidates. Observations at other wavelengths, in particular 
in radio and at very high energies, will most likely be the key to 
clarifying the situation. 
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Appendix A: Uncertainties on gamma-ray 
emissivity estimates for the LMC 

In the estimates of CR relative densities in Sect. 1531 we ne¬ 
glected the difference in metallicity between the LMC and the 
local medium because the corrections due to species heavier than 
He are small (lMorill20()^ . A major uncertainty on the above es¬ 
timates comes from the uncertainty on the gas content of the 
LMC. According to the discussion in lPaner fl. the total gas mass 
in the LMC is quite uncertain because of a possible substantial 
amount of dark gas that is not revealed by the usual radio line 
tracers, but correlates well with the atomic gas dis tribution ac¬ 
cording to infrared observations. Following lPaner J. we adopted 
a total gas mass of (7.2 + 2.4) x 10^ Mq. The uncertainty on the 
gas content may vary from one region to another in the LMC, 
but for simplicity we assumed that all regions have a +33% 
relative uncertainty of their gas content. This directly translates 
into a +33% relative uncertainty of the inferred emissivities. The 
choice of the Xco factor is also formally a source of error, but it 
is minor compared to the uncertainty on the total amount of gas 
(especially since extended emission is found in regions relatively 
poor in molecular gas). 


Appendix B: Energetics and residence time of the 
cosmic-ray popuiation in the LMC 

The emissivities derived in Sect. I5.3l can be used to evaluate the 
total energy of the CR population in the LMC. Emissivity in the 
~ 1 GeV range results predominantly from the decay of neutral 
pions produced by ~10GeV CR protons, and this range of CR 
energies accounts for most of the CR energy density in the local 
ISM. The relative normalisations of emissivities at 1 GeV are 
therefore interpreted as proportional variations of the CR energy 
density in these regions with respect to the local value. 

This estimated energy density is then integrated over volume 
to yield an estimate of a total energy in CRs. The 2D Gaussian 
emissivity profiles were integrated up to 2cr, except for the EO 
component, for which Icr already covers the entire galaxy. The 
physical scale radii were determined from the best-fit angular 
size cr of the 2D Gaussian spatial profile, assuming ?id - 50kpc 
distance to the LMC and correcting for the i - 30° inclination. 
An LMC disk thickness of h - 400pc was used in the calcula¬ 
tion (the s cale height of th e atomic gas was estimated at about 
200pc; see lKim et alJl99^ . Eor the large-scale disk component, 
which has an emissivity three times lower than local, this gives 


Lcr.eo - L4 X 10^^ 


LcR.local 

1 eV cm^^ 


400 pc 


erg, 


(B.l) 


where we assumed a local energy density of ~l-100GeV CRs 
of 1 eV cm“^. This corresponds to the cumulative production of 
CRs by about 1400 supernovae, assuming each explosion re¬ 
leases 10^' erg of kin etic energy, of which about 10 % go into the 


acceleration of CRs (jCaprioli & Spitkovskvll2014 

). The super- 

nova rate in the LMC being about 0.2 SN/century ( 

Hushes et al.l 


l2007h . these 1400 supernovae imply an accumulation of ~1- 
100 GeV CRs over a typical lifetime of 0.7 Myr. This is much 
less than the estimated 10-15 Myr mean res idence time of ~1- 
100 GeV CR protons in the Milky Way disk (iBrunetti & Coding 
I2000h . 

To first order, this factor ~ 10-20 difference on the residence 
time of CRs in the Milky Way and LMC is consistent with CR 
density estimates. The large-scale peak density in the LMC was 
found to be 30% of the local Galactic density, which means 
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about 20% of the peak Galactic density (there is a Galactocentric 
gradient in CR density and t he peak value is ~1.5-2.0 times the 
value at the solar circle; see I Ackermann et^ 1201 1(3) . There is 
therefore a factor ~ 5 difference in large-scale CR density, and 
not 10-20 as might be suggested by the difference in residence 
time. The factor ~2-\ difference can be accounted for from the 
fact that the star formation rate per unit volume is a factor of 
about 3 higher in the LMC than in the Galax\P^. It is also in¬ 
teresting to note that the ~1 Myr residence time for CRs in the 
LMC disk contrasts with the typical age of 5-10 Myr of super¬ 
giant shells, if the latter are confining CRs in localised regions, 
as suggested above. 

Performing the same calculation for the other regions in our 
emission model results in a total CR energy of 1.5 x 10^^ erg, 
3.2 X 10^^ erg, 1.8 x 10^^ erg, and 2.2 x 10^^ erg for components 
El, E2, E3, and E4, respectively (with the same dependencies on 
the local CR energy density and LMC disk thickness as in Eq. 
IB.lb . In total, this makes about 2.3 x 10^^ erg in ~l-100GeV 
CRs, or the equivalent of the production of 2300 supernovae. 
This is probably a lower limit since we neglected the contribu¬ 
tion of additional higher energy particles in these regions fea¬ 
turing a hard spectrum, and the amounts of CRs present at high 
altitude above the disk. 
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Using as star formation rates 2 and 0.2 Mq yr ', and as disk radii 20 
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